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Abstract 

We give an introduction to discrete functional analysis techniques 
for stationary and transient diffusion equations. We show how these 
techniques are used to establish the convergence of various numerical 
schemes without assuming non-physical regularity on the data. For 
simplicity of exposure, we mostly consider linear elliptic equations, and 
we briefly explain how these techniques can be adapted and extended 
to non-linear time-dependent meaningful models (Navier-Stokes equa¬ 
tions, flows in porous media, etc.). These convergence techniques rely 
on discrete Sobolev norms and the translation to the discrete setting 
of functional analysis results. 
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1 Introduction 

A number of real-world problems are modelled by partial differential equa¬ 
tions (PDEs) which involve some form of singularity. For example, oil engi¬ 
neers deal with underground reservoirs made of stacked geological layers with 
different rock properties, which translate into discontinuous data (permeabil¬ 
ity tensor, porosity, etc.) in the corresponding mathematical model. Another 
example from reservoir engineering is the modelling of wellbores; the relative 
scales of the wellbores (~10-20cm in diameter) and the reservoir (~l-2km 
large) justifies representing injection and production terms at the wells by 
Radon measures [34], The mathematical analysis of PDEs involving singular 
data is challenging. The meanings of the terms in the equations have to be 
re-thought; classical derivatives can no longer be used, and weak/distribution 
derivatives and Sobolev spaces must be introduced [4], Beyond these now 
well-known tools, other techniques had to be developed for the most com¬ 
plex models to define appropriate notions of solutions, and to prove their 
existence (and uniqueness, if possible): renormalised solutions [10], entropy 
solutions [3], monotone operators and semi-groups [4], elliptic and parabolic 
capacity [10, 25], etc. The main purpose of this analysis is to ensure that 
the models are well-posed, that is that they make sense from a mathematical 
perspective. It is rarely possible to give explicit forms for, or even detailed 
qualitative behaviour of, the solutions to the extremely complex models in¬ 
volved in field applications. Precise quantitative information that can be used 
for decision-making can be obtained only through numerical approximation. 

The role of mathematics in obtaining accurate approximate solutions to 
PDEs is twofold. First, algorithms have to be designed to compute these 
solutions. But, even based on sound reasoning, in some circumstances algo¬ 
rithms can fail to approximate the expected model [35, Chap. Ill, Sec. 3]. 
Benchmarking (testing the algorithms in well-documented cases) is useful 
to ensure the quality of numerical methods, but it cannot cover all situa¬ 
tions that may occur in field applications. The second role of mathematics 
in the numerical approximation of real-world models is to provide rigorous 
analysis of the properties and convergence of the schemes; this analysis is 
not restricted to particular cases, and is essential to ensure the reliability of 
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numerical methods for PDEs. 

The usual way to prove the convergence of a scheme is to establish error 
estimates; if u is the solution to the pde and Uh is the solution provided by 
the scheme (where h is, for example, the mesh size), then one will try to 
establish a bound of the kind 


|| Uh — w||v: < Ch a (1) 

where || • ||x is an adequate norm and a > 0. Such an inequality provides an 
estimate on the h that must be selected in order to achieve a pre-determined 
accuracy of the approximation. However, major limitations exist: 

• Estimates of the kind (1) can be established only if the uniqueness of 
the solution u to the PDE is known (if (1) holds, then u is unique and, 
actually, the proof of (1) often mimics a proof of uniqueness of u). 

• The constant C usually depends on higher derivatives of u or the pde 
data, and (1) therefore requires some regularity assumptions on the 
solution or data. 

For many non-linear real-world models, including those from reservoir en¬ 
gineering [36] and the famous Navier-Stokes equations, uniqueness of the 
solution is not known unless strong regularity properties on the solution are 
assumed. These properties cannot be established in held applications. Hence 
convergence analysis based on error estimates is doomed to be somewhat dis¬ 
connected from applications. This article presents an introduction to tech¬ 
niques that were recently developed to deal with this issue. These techniques 
enable the convergence analysis of numerical schemes under assumptions that 
are compatible with real-world data and constraints. 

Section 2 details the convergence technique on a simple linear stationary 
diffusion equation. After recalling some basic energy estimates on the model, 
we present the general path (in Section 2.2) to establish the convergence of 
schemes without any regularity assumptions on the data; this path relies on 
compactness techniques and discrete functional analysis tools, translations 
to the discrete setting of functional analysis results pertaining to functions of 
continuous variables. Section 2.3 shows on two particular schemes (two-point 
finite volume scheme, and non-conforming P 1 finite element scheme) how this 
path is applied in practice. In Section 3 we discuss the extension of this con¬ 
vergence technique to non-linear and non-stationary models, more realistic 
representations of physical phenomena. We briefly show that virtually no 
adaptation is required from the technique used in the linear setting to deal 
with the simplest non-linear models. We then give a brief overview of phys¬ 
ical models whose numerical analysis was successfully tackled using discrete 
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functional analysis tools. These include the Navier-Stokes equations, pdes 
involved in glaciology, models of oil recovery, models of melting materials, 
etc. 


2 Convergence by compactness techniques 

2.1 Model and preliminary considerations 

Let us consider, for our initial presentation, the linear diffusion equation 

f — div(AVtt) = / in Q, 

1 u = 0 on dCt. 


In the context of reservoir engineering, (2) corresponds to a steady single¬ 
phase single-component Darcy problem with no gravitational effects [11]; u 
is the pressure and A is the matrix-valued permeability field. This field is 
usually considered piecewise constant (constant in each geological layer), and 
it is therefore discontinuous. Equation (2) cannot be considered under the 
classical sense - with div and V denoting standard derivatives - and must be 
re-written in a weak form; this form is obtained by multiplying the equation 
by a test function v which vanishes on dCl and by using Stokes’ formula [4]: 

{ Find u G Hq(Q) such that: 

VuGilo(D), f A(x)'Vu(x) ■ Vv(x)dx — f f(x)v[x)dx. ^ 
Jn Jn 

Here, Hq(Q) is the Sobolev space of functions v G L 2 (D) (square-integrable 
functions, equipped with the norm ||u||| 2 m) = f Q \v(x)\ 2 dx), that have a 
weak (distribution) gradient Vw in L 2 (f!) d and a zero value (trace) on dCl. 
Under the following assumptions, all terms in (3) are well-defined: 

D is a bounded open set of R d (d > 1) and / G L 2 (Cl), (4) 

A : fl t—>■ Af(i(K) is a measurable matrix-valued mapping, 

30 < a < a < oo such that |H(x)^| < a|£| and H(x)^ • £ > a|£| 2 (5) 

for almost every iGll and all £ G R d . 


Here, j • j is the Euclidean norm on R d . By taking v — u in (3) and by 
applying Cauchy-Schwarz’ inequality on the right-hand side, we find 


®ll |Vm| HUyo) < / A(x)'Vu(x) ■ Vu(x)dx — / f(x)u{x)dx 


< ||./IU 2 (Q)IMlL 2 (n)- (6) 
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Essential to the analysis of elliptic equations is Poincare’s inequality: 

\/v e ifo(ft), |M| L 2 (n) < diam(fi)|| |Vu| Hz^n)- (7) 

Substituted into (6), this inequality leads to the following energy estimate, 
in which the left-hand side defines the norm in //,] (Q): 

||«||Hi(n) : = II |Vu| IU 2 (n) < diam(fi)a“ 1 ||/|| L 2 (n) . (8) 

2.2 General path for the convergence analysis 

Estimate (8) shows that H ( \ (Q) is the natural energy space of Problem (2). 
This estimate is at the core of the theoretical study of (2) and its non-linear 
variants, partly due to Rellich’s compactness theorem [4], 

Theorem 1 (Rellich’s compact embedding) If Cl is a bounded subset of 
R d , d > 1, and if(v n ) ne ^ is bounded in Hq{Q), then (u n )neN has a subsequence 
that converges in L 2 (Cl). Furthermore, any limit in L 2 {fX) of a subsequence 
of (v n ) n£N belongs to H^(Cl). 

This theorem justifies the general path for a convergence analysis that is 
applicable without smoothness assumption on the data or the solution, and 
that can be adapted to non-linear equations. As described by Droniou [14], 
this path comprises three steps: 

1. Establish a priori energy estimates similar to (8) on the solutions to the 
scheme, in a mesh- and scheme-dependent discrete norm that mimics 
the Hq norm, 

2. Prove a compactness result, discrete equivalent of Theorem 1: if (uh)h 
is a sequence of discrete functions that are bounded in the norms in¬ 
troduced in Step 1, then as the mesh size h goes to zero there is a 
subsequence of {uh)h that converges (at least in L 2 (fi)) to a function 
u e Hq(Q), 

3. Prove that if u € Hq (Q) is the limit in L 2 (Cl) as h —> 0 of solutions to 
the scheme, then u satisfies (3). 

Remark 2 The existence of a solution to the PDE does not need to be known. 
It is obtained as a consequence of the convergence proof. 
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The discrete H ( \ (Q) norm is dictated by the scheme. It must be a norm for 
which (i) a priori estimates on the numerical solutions can be obtained, and 
(ii) the compactness result in Convergence Step 2 holds. There is however a 
norm applicable to a number of numerical methods. Let us assume that fl is 
polytopal (polygonal in 2D, polyhedral in 3D, etc.), and that A4 is a mesh of 
D made of polytopal cells. We denote by Hm = madiam(iC) the size 
of A4, and by X M the space of piecewise constant functions in the cells. We 
identify v G X M with the family of its values (vk)kcM i n the cells. £ M is 
the set of all faces of the mesh (edges in 2D), and |er| denotes the (d — 1)- 
dimensional measure of a face cr (i.e. length in 2D, area in 3D). We take one 
point Xk in each cell K, and we let dx,a = dist(xA', cr) (see Figure 1). If o 
is an interface between two cells K and L, then we define d a = dx, a + 
otherwise, d a = d Ki(T with K the unique cell whose o is an face. 



Figure 1: Notations associated with a polytopal mesh. 

A discrete Hq norm on X M is defined by 

imi^-Ei-k(^) 2 - o) 

(tCEm ' / 

Here, and in subsequent similar sums, we use the convention that K and L 
are the cells on each side of o, and that vl — 0 if o C dCl is a face of K. 
This choice accounts for the homogeneous boundary conditions on dQ. 

The major interest of the discrete H q norm, in view of the convergence 
steps 1-3, is apparent in the two following theorems, proved by Eymard et 
al. [30]. Theorem 3 is the key to reproduce at the discrete level the se¬ 
quence of inequalities (6)-(8) leading to the energy estimates mentioned in 
Convergence Step 1 this requires suitable coercivity properties of the scheme. 
Theorem 4 covers Convergence Step 2. Convergence step 3 is more scheme- 
dependent, and relies on consistency and limit-conformity properties of the 
scheme. Theorems 3 and 4 are examples of discrete functional analysis re¬ 
sults. 
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Theorem 3 (Discrete Poincare’s inequality) Let XI be a mesh of 
and set 


9 m = max 


d I\. a 
d'L y cr 


o G Em ) K, L ce ll s on each side of o 


( 10 ) 


If 9 > 9m, then there exists C\ only depending on 9 such that for any v G Xm 
we have \ \v\\ L 2 (n) < Ci\\v\\ H i M . 

Theorem 4 (Discrete Rellich’s theorem) Let (.M n ) ne N be a sequence of 
discretisations of ft such that (9 Mn )neN Is bounded and hM n —t 0 as n —> oo. 
If v n G Xm„ is such that (|]f n ||_ffi,>i n )neN is bounded, then (u n )neN is rela¬ 
tively compact in A 2 (fl). Furthermore, any limit in A 2 (12) of a subsequence 
of(v n ) ne N belongs to 


2.3 Examples 

Besides Theorems 3 and 4, an important feature of the discrete norm (9) 
is its versatility; it is suitable for numerous schemes, even with degrees of 
freedom that are not cell-centred. Here we give a practical illustration, using 
two methods, of the usage of Convergence Steps 1-3 and of the discrete norm 
( 9 ). 


2.3.1 Two-point flux approximation finite volume scheme 

The two-point flux approximation (tpfa) scheme for (2) is given by flux 
balances (obtained by integrating (2) over the cells), and a finite difference 
approximation of the flux — f A(x)Vu(x) -n x(x)dx using the two unknowns 
on each side of cr. 



f(x)dx, 


VA G XI , Vcr G £k '■ Fk,<j — T a {U k — Ul )• 


( 11 ) 

( 12 ) 


Here, £k is the set of faces of a cell K G XI, and the transmissivity r a G 
(0, oo) depends on A and the local mesh geometry [28]. Under usual non¬ 
degeneracy assumptions on the mesh, there exists C 2 > 0 only depending on 
a and a such that 


( 13 ) 
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Convergence Step 1 The inequalities (6)-(8) that lead to the a priori 
estimates on u are obtained by the following sequence of manipulations: (i) 
multiply (2) by v — u and integrate the resulting equation, (ii) apply Stokes’ 
formula, and (iii) use Poincare’s inequality. Since the flux balance (11) is the 
discrete expression of (2), we reproduce these manipulations at the discrete 
level. 

(i) Multiply mid integrate: we multiply (11) by vk = uk and we sum on 
K G M. Accounting for (12) this gives 

E E t <t{ u k ~~ u L )u K = E / f( x )dx u K = / f{x)u{x)dx. (14) 

KcMaCSn K ■' K dfl 

(ii) Apply Stokes ’ formula: this consists of gathering by faces the sum in 
the left-hand side of (15). The contributions of a face are t,j(uk—ul)uk 
and T a [uL — uk)ul = —r a (uK — ul)ul■ Hence, using (13) and Cauchy- 
Schwarz’ inequality on the right-hande side, we find 

C*2 E — ul ) 2 < E t °(. u k ~ u l ) 2 ^ ||/||L 2 (n)|| , w|U 2 (n)- (15) 

(iii) Use Poincare’s inequality: the left-hand side of (15) is C 2 ||w| \ 2 Rl M - 
Invoking the discrete Poincare’s inequality (Theorem 3), we find C 3 
only depending on an upper bound of 9m such that 

< C 3 \\f\\ L 2 (ny (16) 

Estimate (16) is the discrete equivalent of (8) for the solution of the 
TPFA scheme. 

Convergence Step 2 This step is straightforward from (16) by using the 
discrete Rellich’s theorem. This estimate shows that if (A4 n ) ne N is a sequence 
of meshes as in Theorem 4 and if u n is the solution of the tpfa scheme on 
A4 n , then (||^n||ir 1 ,x n )neN remains bounded. Hence, up to a subsequence, 
u n converges in L 2 (fl) towards some function u G Hq(Q). 

Convergence Step 3 As mentioned above, proving that u is the solution 
to (3) hinges on adequate consistency properties enjoyed by the scheme. Here, 
it all comes to the proper choice of transmissivities r a , and to the geometry of 
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the mesh. By taking tp G C£°(f2), multiplying (11) for M. = M. n by p( xk ), 
and summing over all K we hnd 

E E T A( U ^)k - {u n ) L \p(x K ) = E / f( x )<P( x K)dx. 

KeMn <r&e K KeMn J K 

We then gather the sums in the left-hand side by terms involving ( u u )k■ 

E E hr A') V{X L )\ ^ ^ I f {pt2)ip{xK)dx (17) 

KCMn ctCSr KeM n ^ K 

where <p(xl) = 0 if a G S K lies on <9f2. The choice of T a , the geometrical 
assumptions constraining the meshes for the TP FA method (that is, an or¬ 
thogonality requirement of (. XkXl ) and a for a scalar product induced by 
A -1 ), and the smoothness of tp ensure that Ylae£ K t °Vp{ x k) ~ <p{ x l)\ = 
— f K div(AV(p) + \K\0(h,M n ), where \K\ is the d- dimensional measure of K. 
Relation (17) thus gives 


- / u n (x)div(AV(p)(x)dx + 0(\\u n \\ L i(n)hM n ) = / f(x)p(x)dx + 0(h M „ 


in 


where we used the smoothness of <p in the right-hand side. By the convergence 
of u n to u in L 2 (f2), in the limit n —> oo we hnd that u G Hq(Q) satishes the 
following property, classically equivalent to (3): 


V<p G U/°(f2), — f u(x) div (AVip)(x)dx = f f(x)<p(x)dx. 

Jn Jn 

Remark 5 The above reasoning apparently only shows the convergence of 
a subsequence of (u n ) ne pj. However, since there is only one possible limit 
(namely, the unique solution u to (3)), this actually proves that the whole 
sequence (u n ) n eN converges to u. 


2.3.2 Non-conforming P 1 finite element 

Usage of the discrete norm (9) is not limited to numerical methods with 
only/primarily cell unknowns. Let us consider a triangulation T of 2D poly¬ 
gonal domain fl (what follows also generalises to tetrahedral meshes of a 
3D polyhedral domain). The non-conforming Crouzeix-Raviart P 1 finite ele¬ 
ment [9] for (2) has degrees of freedom at the midpoints (x a )ae£ T of the 
triangulation’s edges. The discrete space Yj- of unknowns is made of families 
of reals u = ( u a ) ae g T , where u a = 0 if o C <9f2. These families are identified 
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with functions u : —> M that are piecewise linear on the mesh, with values 

(u a ) a£ £ T at ( x a )ac£ T ■ The non-conforming P 1 approximation of (3) is 


Find u eYt such that: 

Vu G Yp , [ A(x)Wbu(x) ■ Wbv{x)dx 


f(x)v(x)dx 


(18) 


V Jn Jn 

where Vj, is the broken gradient: ( Vbu)\K is the constant gradient of the 
linear function u in the triangle K G T. 


Convergence Step 1 To benefit from Theorems 3 and 4, we need to 
introduce the norm (9), which requires some choice of cell unknowns. Here, 
the most natural choice is to set u as the value of u at the centre of gravity 
x k of K] since u is linear in K. this gives 

WK G T, u K = u(x K ) = - ^ u a . 

ct££k 

This choice associates (in a non-injective way) to each u G Yj- a u = 
(uk)kgt £ Xt- Two simple inequalities, both based on the linearity of 
u inside each triangle, will be useful to conclude Convergence Step 1. 

Lemma 6 Let qr be the maximum over K G T of the ratio of the exterior 
diameter of K over the interior diameter of K. Assume that fj > Then 
there exists C 4 only depending on fj such that, for all u G Yp, 

ll w l|j?i,r ^ C 4 II |Vyu| \\ l 2 ( q ), (19) 

||« ~~ u \\ l 2 (Q) < hr\\ |V&tt| ||l 2 (q>. ( 20 ) 

Proof: Start with (19). There exists C 5 only depending on fj such that for 
all cr G /C we have dist(x^, x a ) < Csd a . Hence, since u is linear inside each 
triangle, 

I % - U L 1 < c \u(x K ) - u(x a ) I c I u(x L ) - u(x a ) I 

d a ~ 5 dist(a;/<, x a ) 5 dist(a;L, x a ) 

< C§ | (V b^)\K | T C 5 1 (Vfell) |a | ( 21 ) 

By squaring (21), multiplying by \o\d a , summing over the edges and us¬ 
ing \ a K — Cq\K\ with Cq only depending on fj, we obtain (19). 

The proof of (20) is even simpler and follows directly from the fact that 

u(x) — u(x) = u{x K ) — u(x) = ( Xbu)\K ' (xk — x) for all x G K. A 
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Equipped with (19) and (20), we now delve into Convergence Step 1. 
Substituting v — u in the formulation (18) of the scheme, the coercivity of 
A entails 

all \^bu\ ||| 2(n) < ||/||l 2 (0)IM|l 2 (Q)- 

Using (20) and hj- < diam(f2), this gives 

a|| |V 6 «| ||| 2(n) < ||/|| L 2 (o)(p||L 2 (n) + diam(U)|| |V 6 «| ||l 2 ( 0 ))- 

A bound on rjj- implies a bound on 6 q- (defined by (10)). Hence, the discrete 
Poincare’s inequality (Theorem 3) and (19) lead to 

II |V 6 u| ||l 2 (0) < {C\C A + cliam(H))a _1 ||/|| L 2 (Q) . (22) 

Estimate (22) is the discrete equivalent of the energy estimate (8). In con¬ 
junction with (19) it gives 

||«|| Hl,T < C A {C 1 C i + diam(n))a _1 ||/|| L 2 (n) . (23) 

Convergence Step 2 This is similar to the same step in the TP FA method. 
If (T n )nen is a sequence of uniformly regular triangulations whose size tends to 
zero, then combining (23) (with T = T n ) and the discrete Rellich’s theorem 
(Theorem 4) shows that u n — > u in L 2 (fl) up to a subsequence, for some 
u G Hq(Q). Moreover, by (20) and (22), we also have u n —$■ u in L 2 (fl). 


Convergence Step 3 Assume now that 

VbU n —> Vu weakly in L 2 (fl) d as n —> oo. (24) 

For p> G C™(Q) we define the interpolant v n G Yj- n by (v n ) a = (p(x a ). The 
smoothness of ip ensures that v n —> tp in L°°(Q) and Vf ,v n —> V(p in L°°(Q) d . 
The convergence (24) therefore allows us to pass to the limit in (18) written 
for u n and v n . We deduce that u G //,] (O) satisfies f Q AVu-Vpdx = j Q f<pdx 
for all smooth ip, which is equivalent to (3). 

The proof of (24) relies on well-established techniques. By (22) the se¬ 
quence (ybU n ) n m is bounded, and therefore converges weakly in L 2 (Q) d to 
some %, up to a subsequence. We just need to prove that x — V?Z. Take 
i/? G C™(Cl) d and, by Stokes’ formula in each triangle, 


VbU n (x) ■ i il>(x)dx = / VbU n (x ) • x/)(x)dx 

KCT n dK 

= £ f K)l K(x)n K . Mx)dS( x ) - W Un(x)divip(x)dx 


KCTn 


I OK 


KeT„ 
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= Z n — / u n (x)divip(x)dx, (25) 

Jn 

where n k is the outer normal to K and (u n )\K denotes values on o from K. 
Since ip — Q on dEl and ip ■ n k + ip ■ n L = 0 on the interface cr between K 
and L, we have 

E E f ( u n)*nK ■ 'ipix)dS(x) 

I<&Tn cr £S K ^ 

= E (u n ) a (n K ■ ip(x)+ n L ■ ip(x))dS(x) = 0. 
ae£,acn 

and thus 

Zn = E E / [(u n )\K(x) ~ (u n ) a \ U K ■ 1 p(x)dS(x). 

KeTn <ree K J a 

By dehnition of {u n ) a we have J a [(u n )\ K (x)-(u n ) a ]dS(x) = 0. Using \{u n )\ K ~ 
(u n ) a \ < diam(/l)|(V6U n )|A'| and the smoothness of ip, we infer 


(K) \k{x) - (u n ) a ]n K ■ [ip(x) - ip(x a )]dS(x) 

< C^ilMn E E W\h K \(V b U n )\ K \ < SC^ClllMn II |VfeU n | || L!(Q) 

with C 7 not depending on n (we used the regularity assumption on T n to 
write \o\hji < C 7 \K\). Invoking the discrete energy estimate (22), we deduce 
that Z n —> 0 and we therefore evaluate the limit of (25) since u n —> u in 
L 2 {Vt) and W h u n —> x weakly in L 2 (Q) d . This gives f n x( x ) • ip(x)dx = 
— f n u(x)divip(x)dx, which proves that x = Vu as required. 



3 Extension to non-linear models 

The previous technique, based on the convergence steps 1-3 and on the dis¬ 
crete Rellich’s theorem and the discrete Poincare’s inequality, would not be 
very useful if it only applied to the linear diffusion equation (2). Conver¬ 
gence of numerical methods for this equation is well-known, and best ob¬ 
tained through error estimates. The power of the compactness techniques 
presented above is that they seamlessly apply to non-linear models, includ¬ 
ing models of physical relevance such as oil recovery and the Navier-Stokes 
equations. Presenting a complete review of these techniques on such models 
is beyond the scope of this article, but we can give an overview of some of 
the latest developments in this area. 
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3.1 Stationary equations 

3.1.1 Academic example 

We first show with an academic example how to apply the previous techniques 
to a non-linear model. We consider 

f -div(A(-,«)V«) = F(u) in Q, ^ 

( u = 0 on dEl ' 

where F : M H > M is continuous and bounded, and A : 12 x M !->• Afd(M) is a 
Caratheodory function (measurable with respect to x e fl, continuous with 
respect to s G M) such that for all set the function A(-, s) satisfies (5) with 
a and a not depending on s. The weak form of (26) consists of (3) with f(x) 
and A(x) replaced with F(u{x )) and A(x,u(x)), respectively. 

As in the linear model case, establishing the convergence of a numerical 
method for (26) by using discrete functional analysis techniques consists of 
mimicking estimates on the continuous equation. Here, these estimates are 
obtained as for the linear model; substituting v — u in the weak form of (26) 
and using the coercivity of A, the bound on F and Poincare’s inequality, it 
is seen that u satisfies 

||«||jfi(n) < diam(n)a _ 1 |n| 1 / 2 ||F|| L oo (K) . 

Writing a numerical method for (26) using a method for the linear equa¬ 
tion ( 2 ) is usually quite straightforward: all f(x) and A{x) appearing in the 
definition of the method (e.g. through r CT for the TPFA method) have to be 
replaced with F(u(x )) and A(x,u(x)), where u is the approximation sought 
through the scheme. A quick inspection of Convergence Steps 1 in Sections 

2.3.1 and 2.3.2 shows that the discrete energy estimates (16), (22) and (23) 
hold with ||/||l 2 (q) replaced with Ifil 1 /-)|F|| L oo (M) . 

Convergence Step 2 then follows from Theorem 4 exactly as in the linear 
case, and we find u G H f ] (O) such that up to a subsequence u n u in L 2 (Q). 
This ensures that F{u n ) —>■ F(u) in L 2 (Q), and that up to a subsequence 
A(-,u n ) —> A(-,u) almost everywhere while remaining uniformly bounded. 
These convergences enable us to evaluate the limit of the scheme by following 
the exact same technique as in Convergence Steps 3 for the linear model. This 
establishes that u is a weak solution of (26). 

Remark 7 Although the strong convergence of u n to u is not necessary in 
the linear case (weak convergence would suffice), it is essential for non-linear 
models such as (26). Indeed, if (it n ) n eN only converges weakly, then F{u n ) 
and A(-,u n ) may not converge to the correct limits F(u) and A(-,u). 
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3.1.2 Physical models 

As mentioned in the introduction, the strength of a convergence analysis via 
compactness techniques is that it applies to fully non-linear models that are 
relevant in a number of applications. 


Elliptic equations with measure data Equations of the form (2) appear 
in models of oil recovery, in which / models wells. The relative scales of the 
reservoir and the wellbores justifies taking a Radon measure for this source 
term [34, 26]. The ensuing analysis is more complex. To start with, the weak 
formulation (3) is no longer suitable [3, 10]. Moreover, due to the singularity 
of the source term, the solution has very weak regularity properties, and 
may not be unique. This prevents any proof of error estimates for numerical 
approximations of these models. 

Discrete functional analysis tools were developed to establish the con¬ 
vergence of the TPFA finite volume scheme for diffusion and (possibly non- 
coercive) convection-diffusion equations with measures as source terms [37, 
24], Key elements to obtaining a priori estimates on the solutions to these 
equations are the Sobolev spaces Wq’ p (EI) (which is Hq(Q) if p = 2), and 
the Sobolev embeddings. The corresponding numerical analysis requires the 
discrete Wq’ p norm on X M 


Wq' v ,M 


£ 


o'I do 


Vk - Vl 

drr 


to generalise the discrete Poincare’s and Rcllich’s theorems to this norm, and 
to establish discrete Sobolev embeddings: if p e (l,d) and q < then 

IMIz^n) < ^\\ v \\wl’ p ,M‘ (27) 


Remark 8 The most efficient proofs of the discrete Poincare’s and Rcllich’s 
theorems actually use the discrete Sobolev embeddings [30, 19]. 


Remark 9 The numerical study of (2) with / measure is currently (mostly) 
limited to the TPFA scheme, since no other method has in general the struc¬ 
ture that enables the mimicking of the continuous estimates [14]. 


Leray—Lions and p-Laplace equations These models are non-linear gene¬ 
ralisations of (2), that appear in models of gaciology [39]. They have a more 
severe non-linearity than (26), since they involve both u and Vw. The gen¬ 
eral form of these equations is obtained by replacing div(AVn) in (2) with 
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div(a(-, u, Vw)), where a : 11 x R x R* 1 4 f 1 satisfies growth, monotony 
and coercivity assumptions. The simplest form is probably the p-Laplace 
equation — div(| Vu| p-2 V«) = / for p G (1, oo). 

Uniqueness may fail for these equations [21, Remark 3.4], which com¬ 
pletely prevents classical error estimates for their numerical approximations. 
Compactness techniques were used to study the convergence of at least three 
different schemes for Leary-Lions equations: the mixed finite volume method 
[13], the discrete duality finite volume method [1], and a cell-centred finite 
volume scheme [29]. These studies make use of discrete scheme-dependent 
Wl' p norms and related discrete Rellich’s and Poincare’s theorems. They 
also require an (easy) adaptation to the discrete setting of Minty’s monotony 
method, to deal with the non-linearity involving Vi. 

3.2 Time-dependent and Navier—Stokes equations 

Studying non-linear time-dependent models requires space-time compact¬ 
ness results. In the context of Sobolev spaces, these results are usually 
variants of the Aubin-Simon theorem [2, 40] which, roughly speaking, en¬ 
sures the compactness in L p (El x (0, T)) of a sequence ( u n ) ne n provided that 
(Vu n ) n£ N is bounded in L p {Ei x (0,T)) d and that (d t u n ) n£ ^ is bounded in 
L q (0,T ; hU _1,r (fl)), where W~ l,r (Vl) = (H / 0 1,r (fl))'. These are natural spaces 
in which solutions to parabolic PDEs can be estimated. 

Carrying out the numerical analysis of these equations with irregular data 
necessitates the development of discrete versions of the Aubin-Simon theo¬ 
rem; this often includes designing a discrete dual norm mimicking the norm 
in W~ 1,r (U). This analysis has been done for various schemes and mod¬ 
els: transient Leray-Lions equations [21], including non-local dependencies 
of a(x,u,Wu) with respect to u (as in image segmentation [33]); a model of 
miscible fluid flows in porous media from oil recovery [6, 7]; Stefan’s model 
of melting material [27]; Richards’ model and multi-phase flows in porous 
media [32]. Discrete Aubin-Simon theorems also sometimes need to be com¬ 
pleted with other compactness results, such as compactness results involving 
sequences of discrete spaces [38], or discrete compensated compactness the¬ 
orems [17] to deal with degenerate parabolic PDEs. 

All these compactness results only provide strong convergence in a space- 
time averaged norm (e.g. L p {Vt x (0, T)) for some p < oo). However, Droniou 
et al. [17, 22] recently developed a technique to establish a uniform-in-time 
convergence result (i.e. in L°°(0, T; L 2 (h2))) by combining the initial aver¬ 
aged convergences, energy estimates from the PDE, and a discontinuous weak 
Ascoli-Arzela theorem. This strong uniform convergence corresponds to the 
needs of end-users, who are usually more interested in the behaviour of the 
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solution at the final time rather than averaged over time. 


Navier—Stokes equations The regularity and uniqueness of the solution 
to Navier-Stokes equations is a famous open problem. Therefore, as ex¬ 
plained in the introduction, the convergence analysis of numerical schemes 
for these equations cannot be based on error estimates. If it is to be rig¬ 
orously carried out under reasonable physical assumptions, this convergence 
analysis can only be done through compactness techniques. 

Let us first consider the continuous case. Because of the term (u-'V)u in 

d\u — A u + (u ■ V)« + Vp = /, (28) 


evaluating the limit from a sequence of approximate solutions (u n )neN re¬ 
quires a strong space-time L 2 compactness on (u n )neN (since {'Vu n ) n£ ^ con¬ 
verges only in L 2 (Q x (0, T))“-weak). Kolmogorov’s theorem ensures this 
strong compactness on (u n ) n eN provided that we can control the space- 
translates and time-translates of the functions. The space translates are 
naturally estimated thanks to the bound on (Vw„) n eN, and the the time 
translates ||w„(- + r, •) — ■UnlUqo, t ; l 2 (Q)) are estimated by 


u, 


i(t+T, x)-u n (t, x)\ 2 dx = 



t+T 


dtu n (s, x)(u n (t + r, x)—u n (t, x))dxds. 


J t 


Equation (28) is then used to substitute d t u n in terms of u n and its space 
derivatives (since divw n = 0, the term involving Vp n disappears). Bounding 
the term ( u n ■ V)w n x u n that appears after this substitution requires Sobolev 
estimates on (u n ) ne pj; these ensure that, only considering the space integral, 
u n G L 6 (Q) and thus |u n |-|Vu n | G L 6 / 5 (h2) (wihout Sobolev estimates, u n G 
L 2 (Q) and |u n | J |Vu n | is not even integrable). 

The same issue arises in the convergence analysis of numerical methods 
for Navier-Stokes equations. Discrete Sobolev estimates of the kind (27) are 
required to estimate the time-translates of the approximate solutions and 
ensure the convergence towards the correct model. Droniou and Eymard 
[16] did this for the mixed finite volume method, and Chenier et al. [8] 
considered an extension of the marker-and-cell (mac) scheme; both refer¬ 
ences establish more scheme-specific Sobolev embeddings than (27), but this 
general inequality is actually sufficient for the analyses carried out in these 
works. 
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4 Conclusions and perspectives 

We presented techniques that enable the convergence analysis of numerical 
schemes for PDEs under assumptions that are compatible with field applica¬ 
tions. In particular, discontinuous coefficients or fully non-linear physically 
relevant models can be handled. These techniques do not require the unique¬ 
ness or regularity of the solutions, and are based on discrete functional anal¬ 
ysis tools - that is the translation to the discrete setting of the functional 
analysis used in the study of the PDEs. 

These discrete tools were adapted to a number of schemes, including the 
hybrid mixed mimetic family [20] (which contains the hybrid finite volumes 
[30], the mimetic finite differences [5], and the mixed finite volumes [15]), the 
discrete duality finite volumes [1], the discontinuous Galerkin methods [12]. 

It might appear from our brief introduction that the discrete Sobolev 
norms and all related results (Poincare, Rellich, etc.) require specific adap¬ 
tations for each scheme or model. This is usually not the case. A framework 
was recently designed, the gradient scheme framework [31, 21, 19], that en¬ 
ables the unified convergence analysis of many different schemes for many 
diffusion PDEs. The idea is to identify a set of five properties that are not 
related to any model, but are intrinsic to the discrete space and operators 
(gradient, etc.) of the numerical methods; convergence proofs of numeri¬ 
cal approximations of many different models can be carried out based on 
these five properties only (sometimes even fewer). Generic discrete func¬ 
tional analysis tools exist to ensure that several well-known schemes - in¬ 
cluding meshlcss methods - satisfy these properties [23], and therefore that 
the aforementioned convergence results apply to these schemes. The gradient 
scheme framework covers several boundary conditions, and also guided the 
design of new schemes [31, 18]. 
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